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ABSTRACT 

It has long been recognised that the main obstacle to accretion of gas onto supermas- 
sive black holes (SMBHs) is large specific angular momentum. It is feared that the gas 
settles in a large scale disc, and that accretion would then proceed too inefficiently to 
explain the masses of the observed SMBHs. 

Here we point out that, while the mean angular momentum in the bulge is very 
likely to be large, the deviations from the mean can also be significant. Indeed, cosmo- 
logical simulations show that velocity and angular momentum fields of gas flows onto 
galaxies are very complex. Furthermore, inside bulges the gas velocity distribution 
can be further randomised by the velocity kicks due to feedback from star formation. 
We perform hydrodynamical simulations of gaseous rotating shells infalling onto an 
SMBH, attempting to quantify the importance of velocity dispersion in the gas at 
relatively large distances from the black hole. We implement this dispersion by means 
of a supersonic turbulent velocity spectrum. We find that, while in the purely rotating 
case the circularisation process leads to efficient mixing of gas with different angular 
momentum, resulting in a low accretion rate, the inclusion of turbulence increases this 
accretion rate by up to several orders of magnitude. We show that this can be under- 
stood based on the notion of "ballistic" accretion, whereby dense filaments, created 
by convergent turbulent flows, travel through the ambient gas largely unaffected by 
hydrodynamical drag. This prevents the efficient gas mixing that was found in the 
simulations without turbulence, and allows a fraction of gas to impact the innermost 
boundary of the simulations directly. Using the ballistic approximation, we derive a 
simple analytical formula that captures the numerical results to within a factor of a 
few. Rescaling our results to astrophysical bulges, we argue that this "ballistic" mode 
of accretion could provide the SMBHs with a sufficient supply of fuel without the need 
to channel the gas via large-scale discs or bars. We therefore argue that star formation 
in bulges can be a strong catalyst for SMBH accretion. 
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1 INTRODUCTION 

The growth of supermassive black holes (SMBHs) at 
the centre of galaxy bulges is known to be correlated 
with observable properties of the host spheroid. The best 
established correlations are between the SMBH mass, Mbh, 
and the velocity dispersion, a, i.e., the Mbh — o" rela - 
tion (|Gebhardt et al. 200CtFerrarese fc Merritt. 2000t ). 

and the Mbh Mbuig e relation 

l|Magorrian et al.. 199dHaring fc Rix. 2004 ). where Mbuigo 
is the mass of the bulge. Understanding the formation and 
growth mechanisms of SMBHs is therefore believed to be 
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important in determining the evolution of the larger host 
systems. 

High luminosities of high redshift q uasars at z ^ 6 sug - 
gest SMBH masses of up to ~ 10^ Mq (|Kurk et al.. 20071 ) . 
implying accretion onto t he black holes close t o or even 
above the Eddington limit (|King fc Pringle. 2007^ for up to 
1 Gyr. Such a strong inflow is most likely the result of major 
mergers between galaxies, which have been shown via hydro- 
dynamical simulati ons to drive gas inward thr ough gravita- 
tional tidal fields (|Barnes fc Hernquist. 199ll ). Simulations 
have also shown that galaxy mergers within large dark mat- 
ter halos at 2 ~ 14 may be able to drive gas in to the cen- 
tre of the galaxy at a suflicient rate to grow a seed BH of 
10* M0 to an SMBH of 10^ M0 in the required time i.e., by 
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z ~ & (|Li et al.. 200?! '). Recent cosmological simulations of 
dark matter halos that include gas physics and black hole 
feedback processes are able to build up the required mass 
by 2 = 6 via exte nded Eddington-limited accretion (e.g., 
ISiiacki et al.. 20091 ). 

However, an ab-initio treatment of all the relevant 
physical processes (gravity, star formation, magneto- 
hydrodynamics, radiative feedback etc.) over the full range 
of scales required is still beyond current numerical capa- 
bilities. Therefore, little is known about the mechanism(s) 
that connect the large-scale flow (~ kpc) to the small-scale 
accretion flow in the immediate vicinity of the SMBH ( ^ 1 
pc). Current cosmological simulations require sub-resolution 
prescriptions to encapsulate the accretion process (e.g., 
ISprineel et al.. 200dBooth fc Schave. 200dSiiacki et al.. 200a ). 
These models commonly use Bondi-Hoyle accretio n 
iHovle fc Lvttleton. 193dBondi fc Hovle. 1944Bondi. 19521 ') . 
corrected upwards by orders of magnitude when the reso- 
lution limitations cause an und er-prediction of the d esired 
SMBH growth rate. Recently, iDeBuhr et al. (20091 ) have 
performed hydrodynamical simulations of major galaxy 
mergers employing a sub-grid model for viscous angular 
momentum transport to dictate the accretion rate, the 
results of which suggest that the Bondi-Hoyle prescription 
actually over-predicts accretion onto the SMBH by up to 
two orders of magnitude. The applicability of the Bondi- 
Hoyle solution to the SMBH accretion rate is therefore not 
clear. 

Looking at the problem from the other end, where 
dynamical and viscous timescales are short compared 
with cosmological timescales, it is expected that ge- 
ometrically thin and op tically thick accretion discs 
ijShakura fc Sunvaev. 197a l form due to a non-negligible 
amount of angular momentum in the gas. While it is 
possible to study such discs at a far higher resolution 
and with more physics, the approach is obviously lim- 
ited by the fact that the disc is commonly assumed 
to just exist a prion. How these discs operate on ~ 
pc scales is also not clear as the discs are suscepti- 
ble to self-gravity and hence may collapse into stars 



cal simulations s how the presence of cold streams (e.g., 
iDekel et al.. 20091 ) immersed into a hotter, lower density 
medium. These complex gaseous structures must also prop- 
agate to smaller scales within the bulge. Secondly, the 
abundance of gas required to fuel the SMBH should also 
fuel star formation in the host. Massive stars deposit large 
amounts of energy an d momentum in the surrounding gas 
(|Leitherer et al.. 19921 ). which is probably one of the domi- 
nant sources of turbulence in the interstellar medi um (e.g., 
iMcKee fc Ostriker. 197iMac Low fc Klessen. 2004 ). In gen- 
eral, where we can observe it, it is evident that tur- 
bulence is ubiquitous in astrophysical flows over a wide 
range of scales and sy stems (for a review see e.g., 
lElmegreen fc Scalo. 2004) . Strong supersonic turbulence is 
a key ingredient in all modern simulations o f star formation 
('e.g.. lKrumholz fc BonneU. 200'?tBate. 20091 ). 

Therefore, we must expect that, far from the SMBH, 
while the specific angular momentum of gas, I, is too large 
to be captured by the SMBH within the central parsec (or 
less), the dispersion in I is also large. It is not unreasonable 
to expect that there will be some gas counter-rotating with 
respect to the mean rotation of the flow. To test the impor- 
tance of these ideas for accretion of gas onto SMBHs, we 
draw on the machinery developed for turbulent flows in the 
star formation fleld by adding a turbulent velocity spectrum 
to the gas in the initial shell. By varying the normalisation of 
the turbulent velocities we essentially control the dispersion 
in the initial gas velocities. 

We flnd that the role of the turbulence and in general 
any disordered supersonic velocity field is two- fold. First of 
all, it broadens the distribution of specific angular momen- 
tum, setting some gas on lower angular momentum orbits . 
Secondly, as is well known (e.g.. iMcKee fc Ostriker. 20071 ). 
turbulence creates convergent flows that compresses gas to 
high densities. We flnd that the dynamics of such regions can 
be described reasonably well by the ballistic approximation 
in which hydrodynamical drag and shocks are of minor im- 
ortance for the high density "bullets" mov ing through the 
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There is thus a genuine gap between these two scales 
and approaches. The purpose of this paper is to try to bridge 
this gap to some degree with the help of numerical hydro- 
dynamical simulations of gas accreting onto a SMBH im- 
mersed in a static, spherical bulge potential. In particular, 
we study the 'intermediate-scale' flow, from the inner 100 
pc of a galactic bulge (just below the resolution limit for 
some of the better-resolved cosmological simulations) to the 
inner parsec, where the gas would be expected to settle into 
a rotationally-supported feature. Instead of assuming the 
presence of the disc, we allow it to form self-consistently 
from the larger scale flow. For simplicity of analysis of the 
results, we start with a spherically-symmetric shell of gas 
of constant density, where the gas is isothermal throughout 
the simulation. The flrst non-trivial element in our study is 
the rotation of the shell, which should force the gas to settle 
into a rotating disc. 

While this turns out to be an interesting system 
in itself, forming a narrow ring rather than a disc, 
we feel motivated to consider a more complex ini- 
tial velocity fleld in the shell. First of all, cosmologi- 



These two effects increase the accretion rate onto 
the SMBH by up to several orders of magnitude, largely 
alleviating the angular momentum barrier problem. 
Note that several authors have already suggested that 
star formation feedbac k and supersonic turbulence in- 
side accretion discs (|Wang et al.. 20091 ) can promote 
SMBH accretion by amplifying angular momentum transfer 
(jCoUin fc Zahn. 200dKawakatu fc Wada. 200dChen et al.. 2009 
An astrophysical conclusion from our simulations is that 
star formation in the host can actually promote AGN 
accretion. 

The paper is arranged as follows. In Section 2 we outline 
the computational method used, and in Section 3 we discuss 
our initial conditions. Section 4 presents an overview of the 
gas dynamics. Sections 5 fc 6 present the results for the no 
turbulence and turbulent cases, respectively, and Section 7 
details the fate of the gas that makes it inside the accretion 
radius. Sections 8 fc 9 comprise our discussion and conclu- 
sion respectively. 
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2 COMPUTATIONAL METHOD 

We use the three-dimensional smoothed particle hydrody- 
namic (SPH)/N-body code G ADGET-3, an updated ver- 
sion of the code presented in ISpringel (2005h . Smoothing 
lengths in the gas are fully adaptive down to a minimum 
smoothing length of 2.8 x 10~^ pc, which is much smaller 
than the scales that we resolve in the simulation. The 
code em ploys the conservative formu lation of SPH as out- 
lined in ISpringel fc Hernguist (20021 ). with the smoothing 
lengths defined to ensure a fixed mass within the smooth- 
ing kernel, for A^'ncigh = 40. Each simulation below starts 
with A^sPH ~ 4 X 10^ particles, with an individual SPH 
particle possessing a fixed mass of meph ~ 12 Mq. The 
Monaghan-Balsara form of artificial visco sity is employed 
l|Monaghan fc Gingold. 198 jBalsara. 19951 ) with a = 1 and 
/3 = 2a. 

The calculations are performed in a static isothermal 
potential with a central constant density core to avoid di- 
vergence in the gravitational force at small radii. We also 
include the black hole as a static Keplerian potential. The 
mass enclosed within radius r in this model is: 

M^f'j — jV/bh + J ' < '"core 

where Mbh = 10* Mq, M„ = 10^" Mq, a = 1 kpc, and 

' core 

= 20 pc and Mcoro = 2 x 10* Mq. The one dimensional 
velocity dispersion of this potential is cr = (GAfa/2a)^''^ = 
147 km s-^ 

As we concentrate on the hydrodynamics of gas accret- 
ing on the SMBH, we make two simplifying assumptions 
which avoid further non-trivial physics (to be explored in 
our future work). First, the position of the SMBH is held 
fixed during the simulations here. This is a reasonable ap- 
proximation as all of our initial conditions are either exactly 
spherically or azimuthally symmetric, or have these symme- 
tries when averaged over the entire simulation volume (e.g., 
turbulence is assumed to be isotropic in this case). Allow- 
ing the black hole to move self-consistently would require 
relaxing the static potential assumption. 

Secondly, gravitational forces between the particles are 
switched off to avoid complications that might arise from gas 
self-gravity. While this could be viewed as a shortcoming of 
our work, we believe that inclusion of self-gravity, ensuing 
star formation and stellar feedback would only strengthen 
our results and conclusions. The gravitational contraction of 
gas clouds would create even higher density contrasts and 
thus make ballistic trajectories even more likely. Star for- 
mation feedback would drive its own turbulence and hence 
amplify the effects of the turbulence that we seed. 

One further simplification is that the gas is isothermal 
throughout the entirety of the simulations. This is a fair 
assumption as cooling times are short for the high densi- 
ties we are considering. Accretion of gas onto the SMBH is 
modelled with the "accretion radius" approach - particles 
that come within an accretion radius of race = 1 pc are re- 
moved from the simulation and we track the total mass and 
each component of the net angular momentum vector of the 
accreted particles. This contrasts with the Bondi-Hoyle ac- 
cretion formulati ons that are frequently used in cosmologica l 
simulations (e.g.. lSiiacki et al.. 200iteooth fc Schave. 20091 ). 
We believe the accretion radius approach, frequently used in 



the star formation field (e.g.. lBate et al.. 19951 ). is essential 
for the problem at hand as it prevents unphysical accre- 
tion of SPH particles with too large an angular momentum, 
whereas the Bondi-Hoyle formulation does not. 

The units of length, mass, time and velocity used in the 
simulations are, respectively, ru = 3.08 x 10 
Mu = 2x 10*^ g (101°Mq), tu = l/n(r^ 
Vn = {GMa/ay/'^ = 208 km s~\ 



cm (1 kpc), 
5 Myr, and 



3 INITIAL CONDITIONS 

The starting condition for all our simulations is that of a uni- 
form density, spherically symmetric gaseous shell centered 
on the black hole. The inner and outer radii of the shell are 
Tin = 0.03 and rout = 0.1 kpc, respectively, for most of the 
simulations (see Table 1 in the Appendix). In principle one 
can expect the outer radius of the shell to be much larger 
in a realistic bulge with effective radius of a few kpc, but 
we are forced to limit the dynamic range of the simulations 
due to computational resources. Furthermore, we believe we 
understand how our results scale with the outer radius of 
the shell (cf. Section l6.2.4p . and hence the dynamic range 
limitation does not actually infiuence our conclusions. 

The total mass of the shell is Msh = 5.1 x 10^ Mq. 
To minimise initial inhomogeneities we cut the shell from a 
relaxed, glass-like configuration. The initial velocity field is 
composed of two parts: net rotation and a seeded turbulent 
spectrum. The net rotation is described by the azimuthal 
velocity component, = Vrot, where Vrot is a parameter 
of the simulation (see Table [T] in the Appendix). In all our 
runs, Drot is below the circular velocity at the shell radii, 
meaning that our initial conditions are not in equilibrium. 
We stress that this is deliberate, as our investigation here is 
the formation of a disc from infalling gas, rather than from 
already rotationally-supported gas. 

Our approach to setting up the turbulent velocity field 
follows that of Dubinski, Narayan fc Phillips (1995). We 
assume a Kolmogorov power spectrum. 



Pvik) 



-11/3 



(2) 



where k is the wavenumber. The key assumption here is that 
the velocity field is homogeneous and incompressible, and so 
we can define v in terms of a vector potential A such that 
V — V X A. This is useful because the components of A 
can then be described by a Gaussian random field with an 
associated power spectrum 



PA{k) . 



(3) 



This is a steep power-law, implying that the variance in 
|yl| will diverge sharply as k decreases, and so we follow 
Dubinski, Narayan fc Phillips (1995) and introduce a small 
scale cut-off fcmin. Equation m can then be written as 



(4) 



where C is a constant that sets the normalisation of the 
velocities. For our purposes we set it equal to unity and 
normalise the velocity field once the statistical realisation 
has been generated. Physically the small scale cut-off fcmin 
can be interpreted as the scale -Rmax — fe"?^^ as the largest 
scale on which the turbulence is likely to be driven. 

Our approach to generating the statistical realisation 
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of the velocity field is straightforward. First we sample the 
vector potential A in Fourier space, drawing the amplitudes 
of the components of Ak at each point {kx,ky,kz) from a 
Rayleigh distribution with a variance given by < \Ak\'^ > 
and assigning phase angles that are uniformly distributed 
between and 2n. Second we take the curl of Ak, 

Vk=ikx Ak (5) 

to obtain the Fourier components of the velocity field. Fi- 
nally we take the Fourier transform of this to obtain the 
velocity field in real space. We use a periodic cubic grid of 
dimension 256'^ when generating the statistical realisation 
of the velocity field and we use tricubic interpolation to es- 
timate the components of the velocity field at the position 
of each SPH particle. 

The initial parameters for each run are given in Table 
[T]in the Appendix. 



4 OVERVIEW OF GAS DYNAMICS AND 
MAIN RESULTS 

Before we embark on a quantitative study of the simulation 
results, we present several snapshots from the simulations 
that illustrate graphically the nature of the gas fiow. In par- 
ticular, here we discuss the overall gas dynamics for two sim- 
ulations that typify the extremes of behavior that we find - 
S30 and S35 (cf. Table [1]). Both simulations have an initial 
rotation velocity Urot = 0.3 which results in a mean circu- 
larisation radius of rcirc = 0.019 (see jjS.lf) . Simulation S30 
has no imposed turbulence. Simulation S35 has turbulence 
characterised by fturb = (I'turbicnt) ~ 1, implying that gas 
turbulent motions are of the order of the velocity dispersion 
in the bulge. 

4.1 A shell with no initial turbulence 

Our no turbulence initial condition acts as a base compar- 
ison for the simulations with seeded turbulence. Figure [1] 
shows the gas column density and the velocity field in the 
angle-slice projection for the simulation S30 at t = 0.06. 
For projection along the z-axis, the gas column density is 
calculated by 

p(x,y)dz , (6) 

where the limits of the integration are given by z{x,y) = 
rtan(", and r = (x^ -f y^Y^"^ . The angle C, is chosen to be 
t&nC, = 1/5 for Figures [1] and [S] This projection method is 
convenient as it permits an unobscured view into the central 
regions where the black hole resides. 

Figure [1] shows the gas flow at an early time, both an 
edge-on (the left hand side panel of the flgure) as well as a 
top-down projection of the shell (the right panel). Due to 
a non-zero angular momentum, gas in the shell is quickly 
pushed aside from the axis of symmetry, opening a cylin- 
drical cavity. Initially gas falls closer to the centre of the 
potential than its circularisation radius, a consequence of 
highly eccentric orbits. As gas moves inside rdrc, centrifugal 
force exceeds gravitational force. Interaction with neighbor- 
ing gas streams results in "radial" shocks in the xy-plane. 



A further set of shocks occurs due to gas initially mov- 
ing supersonically in the vertical direction. As the gas from 
above the z = Q plane collides with the gas falling in from 
below the plane, the particles are shocked and accumulate 
at 2 = due to symmetry, forming a disc. These two sets of 
shocks mix gas with different angular momentum. We shall 
discuss this interesting effect in greater detail below. 

The overall dynamics of the simulation are thus 
relatively simple: angular momentum conservation and 
symmetry dictate the formation of a geometrically thin disc 
in the plane of symmetry of the shell. Accretion of gas onto 
the SMBH would be expected to proceed in an accretion 
disc mode, if at all - realistic models show that star 
formation time scales are much shorter than viscous times, 

depleting the disc of gas before it can feed the SMBH (e.g., 

iGoodman. 200dNavakshin fc Cuadra. 200 jNavakshin et al.. 200 it ). 

We note that our particular choice of velocity field, 
namely a constant azimuthal velocity, when combined with 
a spherical shell, could be viewed as a somewhat peculiar ini- 
tial condition. As we have mentioned, the constant con- 
dition means that the polar regions of the shell are quickly 
evacuated, the gas here spiralling outward in cylindrical ra- 
dius to encounter gas falling in with different angular mo- 
menta, and mixing with the latter. In fact what we are mod- 
elling here is the simplest case of a fiow that undergoes an 
angular momentum re- distribution shock. We shall go into 
more detail on the consequences of this in Section [SI but 
for now we make the point that such a flow is likely to oc- 
cur when gas is accreting in a stochastic fashion from large 
scales, perhaps as the result of a merger. The strong mix- 
ing of the gas with different angular momenta is exactly 
the situation that we wish to explore here, and so we have 
implemented what is essentially a symmetric case of this 
mixing process. We acknowledge that the spherical setup is 
idealised, but it is also the best starting point from which to 
embark on a laboratory of tests where the turbulence and 
rotation is varied. 

4.2 A shell with initial turbulence 

The flow of gas in this case is far more complicated than in 
the case just considered. We flrst show in Figure [5] the full 
slice projection of the simulation volume, e.g., we use here 
x{z,y) — Tout = 0.1 (cf. equation [6]) . The most outstanding 
feature of the flgure are the long dense fllaments that form 
due to convergent turbulent velocity flows. The density con- 
trast between the filaments and surrounding gas is over two 
orders of magnitude by this time. 

Figure [3] shows both the y-axis and z-axis projections of 
this simulation at the same time as [1] did for the no initial 
turbulence case, t = 0.06 . Clearly, some of the filaments 
seen at the earlier time in Figure [2] survive and actually 
penetrate into the innermost region There are now density 
contrasts of as much as three orders of magnitude in regions 
that were completely uniform in simulation S30. We shall 
see that this distinction is a crucial one for the dynamics of 
the gas and SMBH feeding. 

4.3 Turbulence and accretion 

Figure [4] shows the mass accreted by the black hole ver- 
sus time in the simulations S30-S37, e.g., the same rotation 
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Figure 1. Edge-on angle-slice projection (see ^4.1|l of the gas flow at time t = 0.06 in the simulation S30 (left) and top-down projection 
(right). The gas falls in on eccentric orbits, giving rise to a radial shock that propagates outwards as the disc forms. 
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Figure 3. Edge-on angle-slice projection of the gas flow at time t = 0.06 in the simulation S35 (left) and top-down projection (right). 
The velocity field appears far more isotropic that in Figure [T] although an imprint of the imposed net rotation can still be seen. 
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This plot exemplifies the main result of the paper: the accretion 
rate on the black hole strongly increases with increasing levels of 
turbulence when rotation is present (see i|4.3l l. 
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log E[g cm"'] 



Figure 2. Projected gas column density in simulation S35 at time 
t = 0.02, before gas has had a chance to accrete on the SMBH. 
Note the formation of multiple thin and dense filaments due to 
convergent turbulent velocity flows. 



is "chaotic" (|King fc Pringle. 20071 ) rather than large-scale 
disc-dominated. 

We shall spend the rest of the paper investigating these 
results in more depth, suggesting and testing analytical ex- 
planations for the observed trends. 



velocity (and thus angular momentum) but different levels 
of initial turbulence. This demonstrates the main result of 
our study. The accretion rate onto the SMBH strongly cor- 
relates with the strength of the imposed turbulence. The 
accretion rate increases rapidly with increasing Uturb while 
Vtuvh <C Vrot, but then saturates at an approximately con- 
stant level for Vturb ~ frot- The main qualitative explanation 
of the simulations is that turbulence decreases the degree 
to which gas with different angular momentum mixes, and 
creates gas streams with small angular momentum. In par- 
ticular, 

• At low t'turb Vrot, SOI iucrease in the turbulent ve- 
locity leads to greater variations in the density fluctuations 
that are created by the turbulent velocity flows before gas 
circularises (cf. Figure[21). Greater density contrasts decrease 
the amount of angular momentum mixing, resulting in a disc 
rather than a narrow ring. The inner edge of the disc lies 
closer to the accretion radius of the simulations and hence 
feeds the SMBH more efficiently. 

• At I'turb ~ Viot, random initial velocity fields set some 
gas on orbits with a vanishingiy small angular momen- 
tum compared with the mean in the shell. The turbulent 
"kick velocity" in this case almost cancels the mean ro- 
tation for these regions. Since these regions move against 
the mean flow, they are also those that get strongly com- 
pressed. Reaching high densities, they continue to move 
on nearly ballistic trajectories, impacting the innermost re- 
gion on randomly oriented orbits. Accretion in this regime 



5 WITHOUT TURBULENCE: WHY 

"LAMINAR" ACCRETION IS SO LOW 

5.1 Analytical estimates based on circularisation 

of gas 

We shall now argue that the accretion rate of the simulation 
S30, the one with rotation velocity v^ot ~ 0.3 and no initial 
turbulence, is surprisingly low compared with a straightfor- 
ward and seemingly natural theoretical estimate. 

Let us start by estimating the fraction of gas that should 
be accreted by the SMBH in our simulations. Specific an- 
gular momentum of gas determines how close to the SMBH 
it circularises. The angular momentum of a circular orbit is 
given by I = [GM{r)rY^^ . Using equation [T] one can ob- 
tain a general solution for the circularisation radius, rdrc, 
for a given value of specific angular momentum I. We give it 
in two extremes. If rdic ^ ''core, the point-mass Keplerian 
value applies: 



GMbh 



In the opposite 



we have 



GMa 



+ ■ 



4M2 



2Ma 



(7) 



(8) 



which simplifies to rcirc ~ l[a/GMaY^'^ for '"circ 3> r'corc, 
when the second and the last terms on the right hand side 
of the equation are small. 

We now make the simplest possible assumption here by 
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suggesting that gas settles into a disc and that the distri- 
bution of gas in the disc follows the distribution of gas over 
the circularisation radius initially, at time t — 0. Essentially, 
we assume that the Ix and the ly components of the initial 
angular momentum cancel out due to symmetry whereas 
the Iz component is conserved without any exchange with 
neighboring cells. 

To estimate the fraction of gas that will end up inside 
the accretion radius race ^ ^core, we first note it is equation 
[7] that should be used for the circularisation radius of gas 
with a given z-projection of specific angular momentum, l^. 
The requirement h < Ucc = (GMbh^acc)^''^ singles out a 
cylinder with cross sectional radius of ri = y^x^ < 
^acc/urot- The intersection of this cylinder with the shell 
fin < r < rout has volume 27rri (rout — rin). For the shell the 
fraction of the volume that can be accreted is then given 
by the ratio of this volume to the total volume of the shell, 
(47r/3)(rLt-rf„)~(47r/3)rLt: 

r 3G'Mbhracc(rout - Tin) , , 

^■^rot ' out 

For example, for rout = 0.1, rin — 0.03 and Urot ~ 0.1 this 
gives /acc = 0.1, and at t^rot ~ 0.3 we have /acc ~ 10^^ (note 
that G = 1 in code units). 

The latter analytical estimate yields accreted mass 
'7accA/ahcii ~ 5 X 10^ M0, whereas Figure |4] shows that the 
actually measured value to late times is ~ 10^ Mq . The an- 
alytical estimate thus significantly over predicts the amount 
of accretion. 

5.2 Shock mixing of gas: ring formation and end 
to accretion 

Figure [S] explains why our simple analytical theory did not 
work. Here we plot the distribution of gas in simulation S30 
over the ^-component of the angular momentum vector of 
particles, h, at three different times. The initial distribution 
is spread over a broad range of values, with a small but non 
neglibible fraction of gas having Iz < Zacc ~ 0.003, e.g. the 
angular momentum of the circular orbit at rin. In our ana- 
lytical estimate we assumed that this gas accretes onto the 
black hole. Ifowever, this is not what happens. The distribu- 
tion of angular momentum at the later time, t = 0.2, shows 
a strong radial mixing of gas with different angular momen- 
tum. The angular momentum distribution narrows due to 
shocks and eventually becomes a highly peaked Gaussian- 
like ring. 

Note that the deficit of gas in the innermost r ^ 0.01 
region in the second curve in Figure [5] is caused not by the 
SMBH accretion but by the shocks described above. Low 
angular momentum gas shocks and mixes with high angular 
momentum material before it has a chance to travel into the 
SMBH capture region, r < rin. The mixing continues to late 
times and the peak in the angular momentum distribution 
actually moves outwards. This is to be expected as the gas 
that fell in earlier has a smaller angular momentum and is 
in the path of eccentric orbits of gas falling in from greater 
distances. 

We believe this strong mixing of different angular mo- 
mentum orbits is quite a general result of initially "laminar" 
flows. Such flows thus initially form rings rather than discs. 
While our simulations deliberately omit gas self-gravity and 
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Figure 5. The distribution of gas in tiie simulation S30 over the 
2-component of the angular momentum at three different times: 
t = 0.0, t = 0.1, t = 0.2 for solid, long dashed and dashed curves 
respectively. Note how the angular momentum distribution be- 
comes narrower with time (through the action of shocks). 

hence star formation (Sjl]), previous theoretical work shows 
that the fate of the material is then decided by whether 
the viscous time of the ring is shorter than the star for- 
mation consumption time scale. Most authors agree that 
large-s cale gas discs form stars more readily than they ac- 
crete iGoodman. 200jNavakshin et al.. 20071 ). This would 
suggests that "laminar" shells with a finite angular momen- 
tum would form stars more readily than feed the SMBH. 



6 ACCRETION WITH SEEDED 

TURBULENCE: WHY IS IT EFFICIENT? 

6.1 Analytical expectations 

6.1.1 Weak turbulence 

First we consider the case of the turbulent velocity fields 
with nturb ^ Vrot . Figures |S] and [7] show the distribution 
of gas over the angular momentum for three different times 
in the simulations S31 and S32, respectively. As in Figure 
[S] for the no initial turbulence run S30, the initial angular 
momentum distribution is broader than those at the later 
times. The small angular momentum tail of the initial angu- 
lar momentum distribution is more pronounced for S31 and 
especially S32 compared with S30, which helps to explain 
the higher accretion rates measured in these simulations. 

The main effect, however, is the reduction in the shock 
mixing of gas. Consider a "small" angular momentum part 
of the distribution, to be definitive, at Iz ~ 0.01. This part 
of the distribution completely disappears at later times in 
simulation S30, being completely assimilated into the peak 
region. In contrast, in simulations S31 and S32 this part of 
the distribution is only reduced by a factor of 2-3 compared 
with the initial curves. 

Apparently, the initial turbulent velocity flows create 
density contrasts that then propagate to smaller scales. The 
dynamics of a dense region are different from that of a low 
density region, leading to a larger spread in the angular 
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Figure 6. The same as Figure[5]but for simulation S31, with an 
initial turbulent velocity Htm-b = 0.1. There is a greater fraction 
of small angular material in this simulation initially, and more of 
it gets retained in the "tail" to small Iz, e.g., the inner disc, at 
late times. 



(e.g., page 406 in lShapiro fc Teukolskv. 19831 ) for black hole 
accretion. 

For a direct comparison to the simulation results, we 
consider a thin Ar <^ r shell of gas that has angular mo- 
mentum I ~ rViot , and calculate the fraction of gaseous mass 
that has angular momentum small enough to be captured 
within race- In doing so we assume that the energy and an- 
gular momentum of gas are both conserved as the orbits are 
approximately ballistic. The initial specific gas energy, E, 
is small compared with the gravitational binding energy at 
r ~ Tacc, therefore we can set E m 0. An orbit that just 
reaches our inner boundary condition at race has radial ve- 
locity i'r(racc) = 0. This yields the maximum specific angular 
momentum of gas that could still be accreted: 



^max = 2G'Afbhra 



(10) 



The case of a shell that rotates slowly compared with 
the local circular speed, i.e., v^ot ^ = [GMenc ('")/'"] is 
most interesting, since the shell could not rotate faster than 
the circular speed or else the centrifugal forces exceed grav- 
ity; if it rotates at the circular speed then it is rotationally 
supported and would most likely form a disc rather than a 
spherical shell. 
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Figure 7. The same as Figure [6] but for simulation S32. The 
features noted in Figure [6] are even more pronounced here. 



momentum distribution at later times. This increases the 
amount of gas captured by the inner boundary condition. 

6.1.2 Strong turbulence: ballistic accretion 

Now we consider the case when Vturb ~ Vrot. As we have 
seen, the turbulent accretion rate is far larger in our simula- 
tions than that for the "initially laminar" runs. We have also 
seen that differential and chaotic velocity fiows form strong 
density enhancements consistent with the well known result s 
from star formation studies (e.g.. iMcKee fc Ostriker. 2007h . 
High density regions could in principle propagate through 
the mean density gas without much hydrodynamical drag, 
as long as their column densities are much higher than that 
of the surrounding gas. In this case we can approximate gas 
motion as ballistic and use the usual loss cone argument 



Assuming a monochromatic distribution 

As a starting point for this estimate, we approximate the dis- 
tribution of gas velocities in the shell with the turbulent ve- 
locity spectrum by a monochromatic distribution v = Vturb 
randomly and isotropically distributed in directions in the 
frame rotating with rotation velocity i7rot. We emphasise 
here that this is a simple assumption, but a sensible one, 
as on average one would expect the turbulent velocity field 
to have an isotropic character. We relate uturb to the mean 
velocity amplitude in the initial turbulent velocity distribu- 
tion. 

In the simplest case Viot ~ 0, the fraction of solid angle 
that yields angular momentum smaller than Zmax is approxi- 
mately imax/4(r?;turb)^. Therefore, the fraction of mass that 
could end up inside the accretion radius is 

AM GMbh ' acc 



(11) 



where M^h is the mass of the shell. 

However, this approach neglects the rotation of the 
shell. For a non-zero rotation velocity the loss-cone approach 
is valid only for Vturh > i^rot, and equation [TT] must be mod- 
ified to take account of t he orbital m otion. Our derivation 
here is initially similar to lLoeb (20041 ). but taken to second 
order when applying the small-angle approximation. 

We start by considering the axis of zero angular momen- 
tum material for an isotropic, monochromatic distribution 
at a point in orbit of the central black hole. This is defined 
by 



Vvot - "turb COS a = 



(12) 



where a is complementary to the angle between this axis and 
r, as shown in Figure|S] We consider a small perturbation in 
angle about this axis, given by S, which defines the opening 
angle of the loss-cone: 



/ I r-\ ^max 
Viot — Vtuib cos(q ± d) = 



(13) 
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Figure 8. Schematic geometry of an isotropic (in the rest frame of 
the rotating shell) wind, an analogy for the effect of turbulence on 
the loss-cone. The shell rotates with velocity Vrot and the central 
axis of the loss-cone is defined in terms of the angle a. 



We expand this to second-order in the limit that 5 <^ a, 
giving us 

r2 7 

«turb COS a + wturb sm a = (14) 

Z r 

where cos a = Wrot/wturb- We therefore find for the loss-cone 
opening angle the expression: 

1/2 



2lr. 



^'turb sin^ a 



rVrot 



vturh sm a 

Vrot 



(15) 



For an isotropic distribution in the rest frame of the orbiting 
body, the fraction of the solid angle that will be accreted can 
be calculated via 

t:5^ 



flic = 



An 



(16) 



In this case, at a given v^ot, our second-order approximation 
for the loss-cone solid angle is: 

-I 2 



In 



2rvT, 



1/2 



1/2 



(17) 



where ^ = («turb ~ ''rot ) /4i'rot ■ We can express this in the 
limit of two extremes, depending on the relative strength of 
the turbulence with respect to the rotation: 

i) the case when «turb 3> f rot : 

ii) when vtuvh 



/I 7,2 \ 

/ tmax I '-'turb 



+ 



\ 2ri)rot 4i), 

Viot'- 

file 



fturb 



2Vt. 



(18) 



2rVr, 



(19) 



Assuming a Maxwellian distribution 

The treatment we have presented so far is not the whole 
picture, as we have assumed a monochromatic distribution 



for the turbulence. In reality this is not the case. The turbu- 
lent velocity distribution, over all three normally distributed 
components, has a Maxwellian profile, both in the initial 
condition and remaining approximately so in the high den- 
sity regions that form. We find that the distribution for the 
ballistic approximation can therefore be fitted (to within 
~ 10% error) by the expression: 



fiv) 



^?urb^ 



(20) 



To obtain the total mass within the loss-cone we therefore 
need to integrate this expression over the relevant solid an- 
gle: 



Jo 



(21) 



were we should use equation[T7]for the loss-cone solid angle, 
but replacing ^ with ~ {v^ — t^rot ) /4i'rot ■ Unfortunately, 
performing the integration in this case is non-trivial and 
must be done numerically. Luckily, we can obtain a good 
analytical approximation by noting that the modified equa- 
tion [T7] is a maximum when v — Viot, and falls off sharply 
with increasing v. The majority of the accretion will there- 
fore occur in the region of the distribution where v « Urot 
and so we can use equation[l9]in the integral instead, giving 
us: 



Ml, 



2ln 



-(I'/'-turb)' 



dv 



(22) 



where our limits are a perturbation either side of Urot, the 
size of which we choose to be At; = i;rot/4Q. Performing the 
integral yields 



Ml, 



In 



2rVrot 



Sln 



4ri,turb7ri/2 



(23) 



where 



tp = erf 



5«r, 



4i;tur 



erf 



V 4i;turb / 



For a thick shell we should generalise equation [23] by using 
the density weighted average of r over the shell, (r). For an 
initially constant density profile used in our simulations and 
rout ~ 3rin, we have (r) ~ 2rout/3, and hence 

AM _ 3(GMbhracc)'''' ^ 9(GMbhracc)'/' ^^4) 

Mahell 2^rout?;rot 4^/2r^;turb7^1/2 

As before, we can analyse this equation in various limits: 
i) when uturb >• "rot: 

AM _ (2GMbhracc)''" 



Mshell rout?;turb 

u) when VtuA ~ Vrot- 

AM (GMbhra 



1 + ' 



vl/2 



M,h 



(25) 



(26) 



1 it should be noted that this choice is somewhat arbitrary but 
our derivation here is of course approximate 
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iii) when nturb «C Vrot'- 

AM _ 9(GMbhra 



Nl/2 



(27) 



Mshell 4%/2rout?;turb7ri/2 

where we have dropped numerical factors of approximately 
unity. 



^turb = 0, by about 5 orders of magnitude. The role of tur- 
bulence in simulations S50 - S57, as suggested in i)4.3l is 
to increase the accretion rate. This occurs by enhancement 
of the angular momentum distribution into the low angu- 
lar momentum end, and also by weakening the efficiency of 
shock mixing of the gas. 



6.2 Detailed analysis of results 

Our ballistic approximation to the accretion rate yielded a 
result that depends on several parameters of the simulations, 
e.g., the "accretion" radius race, the turbulent velocity tJturb, 
the rotation velocity ii^ot and the dimensions of the shell. We 
shall now look at simulations covering a range in the space 
of these parameters to discuss specific trends and determine 
whether the simulations support our simple analytical the- 
ory. 



6.2.1 No net angular momentum case 

Figure [9] shows the accreted mass (left panel) and the accre- 
tion rate (right panel) as functions of time for simulations 
SOO - S07 (see Table 1). The rotation velocity for these runs 
is set to zero, and the turbulent velocity parameter vt is 
varied from (solid curve) to 2 (blue short dashed). 

The situation is rather simple for these runs. For turbu- 
lent velocity parameters much smaller than unity, the flow 
hardly deviates from that of a freely falling shell with a neg- 
ligible pressure. Due to the constant density profile in our 
initial shell, this naturally leads to the accretion rate in- 
creasing as M oc t^ (see the dotted power law in the right 
panel of Figure[9]) starting from the free fall time of the inner 
shell to that of the outer shell. As turbulence increases to 
t'tuib ~ 1, however, two effects are noteworthy. Firstly, the 
accretion episode starts earlier and also finishes later. This is 
naturally due to the spread in the initial radial velocities of 
the turbulent gas, with some regions starting to fall in with 
negative velocities, and hence arriving at the SMBH earlier, 
and others doing the opposite. There is also a significant 
reduction in the accretion rate at the highest values of tur- 
bulence. We believe this is due to the same effects already 
discussed in M.3\ Convergent turbulent flows create high 
density regions with a large range of angular momentum. 
These regions do not mix as readily as mean density regions 
and thus retain their angular momentum for longer. This re- 
duces the amount of gas accreted. However we note that in 
this rather unrealistic set up (zero net angular momentum) 
gas accretion is still very efficient, with ~ 20 — 100% of the 
shell being accreted by time t = 0.2. 



6.2.3 Rotation and turbulence parameter space 

We now combine the results of all our numerical experiments 
from SOO to S57 in Table 1 by considering a single character- 
istic - the accreted gas mass at time t = 0.2 in code units. 
All of these runs have the same initial geometrical arrange- 
ment of gas and a fixed accretion radius race = 10~^ but 
differ in the strengths of the initial turbulent and rotation 
velocities. 

The results are displayed in two ways. Figure [Til shows 
the accreted mass versus rotation velocity, «rot • The simula- 
tions with the same level of mean turbulent velocity, t^turb, 
are connected by the different lines. In particular, the solid 
black curve shows the results for t'turb = 0, and the dashed 
blue curve shows the results for Wturb = 2. 

It is clear that increased rotation velocity, and thus net 
angular momentum, always reduces the accreted mass, for 
any value of the turbulence parameter. This is in a general 
agreement with analytical expectations explored in i)6.1.2l 

Together with the simulation results, we plot the pre- 
dicted maximum accreted mass fequation l26p . which shows 
excellent agreement with the maximum delineated by the 
levels of turbulence at saturation. The analytical fit here is 
described specifically by 



AM 

Mshcll 



(Gil/bhra 



(28) 



since the maximum fraction that can be accreted is unity 
and for Vrot = we would expect the entire shell to have 
accreted by late times. 

The increased level of turbulence acts to decrease the 
accretion rate in the case of very small rotation velocity 
lirot ~ but increase the accretion rate at higher rotation 
velocities (albeit up to saturation). This is not contradic- 
tory at all, however - the effect of turbulence in all cases 
is to spread the initial angular momentum distribution and 
prevent it from rapidly mixing into a single peak. In the 
case of no initial rotation this reduces the accretion rate by 
moving gas from zero to finite angular momentum orbits, 
whereas in the case of "large" initial net angular momen- 
tum the accretion rate is increased as some gas is moved to 
the low angular momentum orbits. 



6.2.2 The highest rotation velocity case 

In the other extreme. Figure \T0\ shows the accreted mass 
(left panel) and the accretion rate (right panel) for the sim- 
ulations S50 - S57, where the rotation velocity is set to the 
maximum explored in the paper, e.g., Urot ~ 0.5. The range 
of the turbulent velocity parameter and the respective curve 
coding are the same as in Figure [9] 

Comparing Figures [10] and [9l it is obvious that rotation 
significantly decreases the amount of gas accreted, and the 
accretion rate. The reduction is most severe in the case of 



6.2.4 Accretion and the size of the shell 

Equation [21] predicts that the fraction of accreted mass 
should be inversely proportional to the shell's size, which we 
characterise by rout . To test this prediction we repeated the 
simulation S35 for a scaled-up version, rout ~ 0.2, rin = 0.06 
and a scaled-down, rout ~ 0.05, rin = 0.015 shell (see runs 
S60 and S61 in Table 1). 

To be definitive, we compare the accreted gas mass at 
the dynamical time at the outer radius of the shell, which 
corresponds to just after the initial peak in the accretion 
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time [5 Myr] time [5 Myr] 

Figure 9. Accreted mass vs. time (left) and accretion rate vs. time (right) for a velocity field that contains no net rotation (iJrot = 0). The 
different lines represent different strengths of the mean turbulent velocity, ■Wturb- (black solid), 0.1 (black dotted), 0.2 (black dashed), 
0.3 (black dot-dashed), 0.5 (brown dot-dot-dash), 1.0 (red dashed), 1.5 (pink dotted), and 2.0 (blue dashed). Analytical line (thick red 
dashed) displays aMont^ slope. 




time [5 Myr] tl^ig [5 Myr] 

Figure 10. Accreted mass vs. time (left) and accretion rate vs. time (right) for a velocity field that contains high net rotation (?^rot = 0.5). 
Linestyles are as per Figure |9] 



rate. This is approximately the time at which we would 
expect the ballistic mode to come to an end and a disc- 
dominated accretion mode to begin. The results are plotted 
in Figure [12] (left panel) as a function of rout- 

6.2.5 Trend with accretion radius 

The accretion radius is a "nuisance" parameter of our study 
in the sense that it is introduced only to allow the simu- 
lations to run in a reasonable time. For physical reasons it 
would be desirable to make race as small as possible. There- 
fore it is important to test whether the suggested analytical 
scaling of the results indeed holds. With this in mind, we 
repeated same simulation S35 for 5 different values of race, 
spanning a range in race from 0.001 to 0.005. The accreted 
mass from these simulations is shown in Figure [12] (right 
panel). The red dotted power-law gives the analytical pre- 



diction. Although similar, it is somewhat less steep than the 
simulation results. 

6.2.6 Visualisation of the circumnuclear disc 

Figure [T51 compares the surface density projected along the 
z-axis for simulations S30 (left) and S35 (right) in the inner 
r = 0.04 region. Figure [TJ shows the same for these two 
snapshots but projected along the y-axis. 

As was already clear from Figure[5l the gas forms a nar- 
row ring in the simulation S30. This occurs due to efficient 
mixing of material with different angular momentum. This 
mixing is both radial (the ring is narrower) and vertical (the 
ring is vertically less extended). 

In contrast, in the simulation S35, turbulence creates 
a broad distribution of angular momentum which does not 
mix so well. Mixing in the vertical direction is actually rea- 
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Figure 11. The accreted mass trend with rotation velocity, plotted for different levels of turbulence at t = 0.2 (left) and the accreted 
mass trend with mean turbulent velocity, plotted for different levels of rotation at t = 0.2 (right). 
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Figure 12. Scaling of accreted mass with rout (left) and race (right) for simulations with Vrot = 0.3 and iijurb = 1-0. Symbols correspond 
to the dynamical time at the outer radius of the shell i.e., t = 0.05, t = 0.1, t = 0.2 (left) and t = 0.1 (right). Error bars are somewhat 
arbitrary, as it is not clear when the main ballistic mode of accretion ends - here they denote ±t/3. Normalised analytical fits are equation 
1241 (red), the actual fit that would be expected for the simulation parameters, and equation 1261 (green), the maximum accretion expected 
at saturation. 



sonably efficient, as the disc is thin and lies close to the xy 
plane (see the right panel in Figure [14)) . However the radial 
mixing is not so strong and the disc radial structure is very 
extended in S35 compared with S30. 

6.2.7 Radial structure of the disc 

Our analytical theory also makes a prediction for the sur- 
face density of the disc in the innermost part. Indeed, the 
mass captured by the inner boundary of radius race scales 

1/2 

as Mcapt oc race • This mass would likely sit in a rotationally 
supported disc of roughly r ~ race size. The surface density 
in the disc should then behave as E oc Mcapt (T')/r^ oc r^^^^. 

This analytical prediction (red power-law) is compared 
to the simulation results in Figure [15] which shows the sur- 
face density defined on radial shells for all the different lev- 



els of turbulence for v^ot = 0.5. While the region beyound 
r ~ 0.01 shows separate rings and highly variable surface 
density profiles for different simulations, the disc inside the 
region r ~ 0.01 has a similar radial surface density profile 
for all simulations with uturb ~ 0.5. The agreement between 
simulations and theory is quite reasonable in that region ex- 
cept for runs with a small turbulent velocity t;turb < Vrot. 
For these simulations the ballistic approximation is not ap- 
propriate and the gas does not penetrate into the r < 0.01 
region. At radii larger than 0.01 the surface density distri- 
bution is dominated by the interactions with the bulk of the 
initial shell as it continues to fall in. The ballistic approxima- 
tion is of a limited use here because gas settling in the disc to 
late times (i) went through many collisions/interactions and 
hence the orbits are not ballistic; and (ii) has a relatively 
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Figure 14. Slices of the CND (using tanf = 2/5) formed with no turbulence and with mean Wturb = l-Oi both for tirot = 0.3 at t=0.418 



large angular momentum, so the small angle approximation 
we made in i]5.1lis violated. 



The observed E oc r~^^'^ dependence in the inner disc 
is another piece of simple but convincing evidence that the 
ballistic accretion approximation that we made in i]6.1.2l is 
a reasonable one for the innermost region of the flow. 



6.2.8 Vertical structure 

Figure [14] showed edge-on projections of the disc for the 
simulation S30 (no turbulence) and (S35) initial turbulent 
velocity Wturb = 1. It is clear that the vertical stucture of the 
disc is actually rather complex, especially beyond r ~ 0.02 
region. This is intuitively to be expected as the dynamical 
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Figure 15. Azimuthally-averaged surface density versus radius 
for runs S50-S57 at time t = 0.418. Coloured curves correspond 
to simulations with different levels of mean turbulent velocity, as 
in Figure [9] The red line above the curves shows the theoretically 
expected (un-normaliscd) power law S oc r~^^'^. 



Figure 16. Dependence of vertical disc scalehcight, H/r, versus 
radius, for all levels of turbulence, with v^ot = 0.5 (S50-S57), at 
time t=0.418. Line styles of the curves are the same as in Figure 
[9] The red line above the curves shows the theoretically expected 
result, for the (un-normalised) slope H/r oc r. 



time is longest at larger radii, and therefore circularisation 
and decay of vertical velocity dispersion in the disc should 
take longer here. We can further quantify this by assuming 
that the vertical velocity dispersion is pumped at the rate at 
which the matter is brought in. Generation of momentum 
due to turbulence will therefore be described by pturb ~ 
Eucirc- These motions should decay on the local dynamical 
time, r/vcirc- Assuming a steady-state allows us to estimate 
the vertical velocity dispersion cLS Cy vert ~ Uoutr/rout, where 
Wout is the circular velocity at rout- Therefore, the expected 
disc aspect ratio in the outer parts (where the potential is 
isothermal) is H / R ^ r/rout oc r. 

We can calculate a dynamical, azimuthally-averaged 
disc scaleheight for the circumnuclear disc that forms by 
using the projection of the velocity dispersion along the 'ver- 
tical' direction for each annulus i.e., the mean angular mo- 
mentum vector in each case. We therefore set. 



GMe„c(r) 



(29) 



where Monc(r) is the enclosed mass at r, avert (r) is the 'ver- 
tical' velocity dispersion at r. Figure [TB] plots ratios H/r as 
functions of radius for simulations S50-S57. The solid line is 
for simulation S50, for which no initial turbulence is present. 
It is notable that the disc formed by the flow in this case 
is geometrically much thinner than in all the cases S51-S57. 
Among the turbulent cases, the general trend in H/r is ex- 
actly as predicted above. The large scatter is explained by a 
strong time dependence: the disc is being stirred by strong 
waves that would eventually circularise the disc. 

In the region inside r — 0.01, however, the observed 
dependence H/r is flat while the analytical expectation is 
H/r oc r^/^. We note that the vertical disc scale height is 
much larger than the SPH smoothing length, implying that 
the result inside r = 0.01 is not due to under resolving of 
that region. We believe that the larger than expected vertical 
scale height is due to dynamical heating of the inner disc by 
the waves already noted at larger radii. Indeed, the inner 



disc is far less massive than the outer one, and any waves 
present in the outer region have to either dissipate or get 
reflected in this region. 



7 ACCRETION INSIDE THE CAPTURE 
RADIUS 

While we cannot resolve the hydrodynamics of gas within 
the accretion radius, race (1 pc for most of the simulations 
presented here) , we track both the mass and the angular mo- 
mentum of the particles that are captured. We can therefore 
calculate the mean angular momentum of the gas that set- 
tled inside race. If the disc viscous time is comparable or 
longer than th e duration of our s imulations, which is likely 
fe.g.. iNavakshin" fc Cuadra. 20051 '). then the direction of the 
mean angular momentum of the gas captured within race 
region is related to that of the sub-parsec disc. Note that in 
principle the disc can be warped and then its structure is 
more complex (Lodato et al 2009). 

Figure [T7] shows the evolution of the direction of the 
mean specific angular momentum for the gas particles ac- 
creted inside race in the simulation S30 (left panel, no initial 
turbulence) and S37 (right panel, fturb = 2). The latter sim- 
ulation is initialised with the highest levels of turbulence we 
considered in this paper, and therefore the effects of disc 
plane rotation are the strongest. The "north pole" in the di- 
agram corresponds to a specific angular momentum vector 
oriented along the z-axis. Clearly, there is negligible devia- 
tion from that direction for the "laminar" simulation S30, as 
expected. For the simulation S37, however, there is a com- 
plicated time evolution of the innermost disc direction. At 
one point the disc appears to be completely counter-aligned 
to the mean angular momentum of the shell. 

The accretion disc "plane wandering" at r ^ race should 
also be seen in the larger scale flow of gas. Figure [TS] shows 
the innermost r = 0.005 region for the simulation S37 at 
time t = 0.072. We can clearly see here that the disc in 
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Figure 17. Plots of the accreted angular momentum direction in the Aitoff projection. Both plots are with i^iot = 0.3, with no turbulence 
(left) and Vturb = 2.0 (right). Crosses mark every 0.001 in code time units. Note that the title of "SMBH spin direction" refers to the 
spin orientation that the SMBH would have if it had accreted all of the gas within the capture radius; in reality this is of course not the 
case. 
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Figure 18. Projected gas density and velocity field for simulation 
S37, the innermost r = 0.005 region, for t = 0.072. Note how 
strongly the disc is misaligned with the xy-plane. At later times 
the misalignment decreases as more material accumulates into the 
disc. 



the immediate vicinity of race is also strongly displaced in 
terms of direction from the nominally expected one. One 
should also note the flows of gas along directions roughly 
orthogonal to the disc plane. These flows must have very 
different orientations of angular momentum, and it is these 
that cause the plane of the "sub-grid disc" to execute the 
non-trivial path seen in Figure [TTl 



8 DISCUSSION 

8.1 Summary of the results 

We have studied the hydrodynamics of an initially uniform 
rotating gaseous shell, seeded with an initial turbulent ve- 
locity spectrum as it falls into the inner part of a spherically 
symmetric static potential. Particular attention was paid to 
the amount of gas that made it inside the "accretion radius" 

^acc • 

In the control runs without initial turbulence, gas settles 
into a nearly circular rotationally supported ring (see Fig- 
ures O [131 [Mj. This was somewhat of a surprise. Intuitively 
we expected a wide (Ai? ^ R) disc, as the initial distribu- 
tion of angular momentum is broad (Figure[5]). Furthermore, 
only a few hundred solar masses of gas were accreted, com- 
pared with the ~ 5 X 10* M0 expected on the basis of a 
simple circularisation radius argument (Sj5TT]). In hindsight, 
both of these effects are caused by the effective mixing of 
gas by the shocks. Low angular momentum gas gets mixed 
with high angular momentum gas before the former reaches 
the capture radius r = race- The inner edge of the ring is a 
significant distance from the accretion radius. In the context 
of AGN accretion, such a ring would probably be strongly 
gravitational ly unstable and form stars rather than feed the 
SMBH (e.g.. lGoodman. 200dNavakshin et al.. 2007l \ 

In contrast, runs in which the initial shell was seeded 
with turbulence show qualitatively and quantitatively differ- 
ent results. The most significant of these is that the accretion 
rates are higher by up to 3-4 orders of magnitude compared 
with the simulations that had no initial turbulence imposed. 
This is particularly striking since the net angular momen- 
tum of these simulations is the same. 

It is obvious that turbulence broadens the angular mo- 
mentum distribution, setting some gas on low angular mo- 
mentum orbits. A less trivial conclusion from our numerical 
results is that turbulence also prevents the efficient gas mix- 
ing pointed out above. Physically, high density regions are 
created by supersonic turbulent convergent flows. These re- 
gions can then travel nearly ballistically through the average 
density gas. Some of this gas is able to reach the accretion 
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radius, increasing the SMBH accretion rate compared with 
runs without turbulence. This parameter space trend lasts 
until the initial mean turbulent velocity becomes greater 
than the net rotation velocity, at which point the accretion 
trend saturates and begins to slowly decrease as turbulence 
is increased further. 

The development of high density regions in su- 
personic turbulence is by now a classical result 
from a number of papers (for reviews see e.g., 
lElmegreen fc Scalo. 2004McKee fc Ostriker. 20071 ). There 
is also a body of work that highlights the complicated 
dynamics of sho cks in the multi-phase interstella r medium. 
For example, iBonneU et al. (200dDobbs (20081 ) suggest 
that shocks in a multi-phase medium generate (or rather 
preserve) differential motions, allowing the high density 
regions to slip through the shocks. These effects are similar 
in nature to the ballistic behavior of high density clumps 
fou nd in our simulations. 

ISchartmann et al. (20091 ) simulated young nuclear star 
clusters injecting the energy into the surrounding gas via su- 
pernovae and winds. As in our study it was found that high 
density streams form. Those with small angular momentum 
feed th e small scale "torus" . These results, and earlier re- 
sults bv lCuadra et al. (20061 ) for stellar winds in the central 
parsec of the Milky Way also support the viability of the 
"ballistic" AGN feeding mode. 

We have attempted to build a simple analytical theory 
that would predict the accretion rate for our capture in- 
ner boundary condition. We first calculated the accretion 
rate in the ballistic approximation assuming that gas in 
the shell receives monochromatic velocity kicks v = titurb in 
random directions in the frame rotating with velocity Viot- 
This "monochromatic" theory predicts that gas capture rate 
(fraction) is the highest at Wturb ~ ^'rot. We then took into 
account in an approximate way the corrections associated 
with the fact that the real turbulent velocity distribution 
function is not monochromatic. The fraction of gas accreted 
within race calculated in this way turned out to describe our 
numerical results well. Furthermore, the theory predicted 
that the inner disc surface density should have a power-law 
scaling: E(r) oc r^'^^^, which is indeed born out by the sim- 
ulations. 

The approximate agreement of theory and simulations 
in this paper implies that we can use our results for sys- 
tems of astrophysical interest by an appropriate rescaling of 
parameters, which we attempt to do in the next section. 



8.2 Can quasars and AGN be fed by ballistic 
accretion? 

The most serious challenge to theories of quasar and 
AGN feeding is gravitational instability of these discs 
and the ensuing fragmentation of discs into stars (e.g.. 



phe" of AGN discs can be avoided is to depo sit the gas at 
small enough radii (e.g.. lKing fc Pringle. 20071 ). Namely, one 
can show that gaseous discs are hot enough to avoid grav- 
itational fragmentation within the self-gravity radius, _Rsg, 
which is of the order of 0.01—0.1 pc depen ding on the SMBH 
mass and some of the disc parameters (see lGoodman. 20031 ). 
Inclusion of feedback from star formation allows the disc to 
survive and yet maintain an interesting accretion rate on the 
SMBH out to a few times larger radii (Pedes and Nayakshin, 
in preparation). 

We shall now try to estimate the plausibility of such a 
small radius gas deposition in the context of the "ballistic" 
gas deposition model we have developed. To this end we 
shall use the approximate scalings for the fraction of gas 
accreted from a rotating turbulent shell derived in this paper 
in i mX^ 

The size of the shell, rout, enters this estimate di- 
rectly. Existence of the Mbh — Mbuigo relation implies 
(e.g., iHaring fc Rix. 20041 ) that the SMBH must be fed 
by gas locate d roughly the size of the bulge away from 
the SMBH jKing. 200jKing. 20051 ). or else it is diffi- 
cult to see why the qu oted feedback link must exist 
iNavakshin fc Power. 2003 ). We estimate that the effective 
radius of the bulge, Rb, is about 2.5 kpc for veloc ity dis- 
persion a — 150 km/sec (see iNavakshin et al.. 2003 ). If the 
rotation velocity at Rb is a small fraction of a, say v^ot ~ 50 
km s~^, then the fraction of the gas that can be accreted 
from a turbulent shell is 



O.OOIGA/, 



1/2 



/ 2.5kpc \ / 50km s"^ 
0.2pc J \ rout / V ^rot 



(30) 

where Mg = Mbh/ 10* Mq. Here we have used equation 1261 
assuming that the mean turbulent velocity is of the order 
of the rotation velocity, which is our maximum in the ac- 
cretion trend (although it should be noted that going to 
stronger turbulence does not decrease the accreted mass by 
much - cf. ii6.2|) . The Afbh — Afbuigc relation here gives the 
stellar mass in the bulge as ~ 500Mbh, and if we assume 
that ~ a half of the bu lge gas was used up in star formation 
(|McGaugh et al.. 20091 ). we find that the action of super- 
nova feedback in the context of our model yields an accreted 
mass of ~ 5 X lO'' Mq . This accretion would take place over 
roughly the dynamical time at the outer edge of the bulge, 
giving us an accretion rate ~ a few MQyr"^. In this simple 
estimate, then, the SMBH is indeed able to grow as mas- 
sive as the observations require, assuming that our ballistic 
accretion mode can be sustained for the required time. Fu- 
ture work should test these ideas more directly, with "tur- 
bulence" driven directly by feedback from star formation in 
the bulge, where the bulge gas distribution has been derived 
from initial conditions in the gas on larger scales, e.g., the 
virial radius of the host dark matter halo. 



iPaczvhski. 197dKolvkhalov fc Sunvaev. 198ClCollin fc Zahn. 199dGoodman. 200 jNavakshin et al.. 20071 ). 



Observations (e.g., Paumard et al.. 200dBartko et al.. 20091 ) 

and numerical simulations 

(|Bonnell fc Rice. 200dHobbs fc Navakshin. 20091 ) suggests 
that this process operated in the central 0.5 pc of the Milky 
Way, creating young stars and probably leaving the resident 
SMBH - Sgr A* - without gaseous fuel to shine now (e.g., 
iLevin fc Beloborodov. 200 jNavakshin fc Quadra. 20051 ). 

One obvious way in which the self-gravity "catastro- 



8.3 A positive feedback link of star formation to 
AGN activity? 



Observations show that starburst and 
tivity are correlated in a number of 



AGN ac- 
ways (e.g.. 



ICid Fernandes et al.. 200]lGonzalez Delgado et al.. 200]|Farrah et al.. 20031 ). 
One might think that this is entirely natural, as both phe- 
nomena require a source of gas to be present. However, star 
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formation on kpc-scales does not have to be casually con- 
nected with SMBH activity. As an example of this, consider 
a case when the angular momentum of the gas is large, so 
that a kpc-scale disc is formed. If star formation in this disc 
simply consumes the gas then there is nothing left to feed 
the SMBH. This is what may have happened on smaller 
scales in the central parsec of our Galaxy ~ 6 million year s 
ago jNavakshin fc Quadra. 200dNavakshin et al.. 200?! ). 
Naturally, if such a process can rob the SMBH of fuel at 
parsec scales then the situation becomes even worse at 
kiloparsec scales. 

The opposite situation is also possible because the 
SMBH mass is a small fraction of the bulge mass. It is not 
implausible to have a sub-pc scale, non self-gravitating ac- 
cretion disc that could feed a moderately bright AGN for a 
long (10^ — 10^ years, say) time. The amount of gas involved 
in this could be 10^ — 10^ Mq, i.e., miniscule by galaxy stan- 
dards. Why such activity should always be connected with 
a powerful starburst is not clear. 

However, the inclusion of star formation feedback could 
alter this picture. For example, within a kpc-scale disc, 
feedback from stellar winds and supernovae gives velocity 
kicks to the surrounding gas. For a standard IMF, the to- 
tal momentum in supernova shells and winds from mas- 
sive stars, integrated over their life time, is about e*Af*c, 
where e « 10~^, and M, is the total mass of the ste llar 
population (jLeitherer et al.. 199!frhompson et al.. 20051 ). If 
this momentum is absorbed by gas with mass Mg then the 
average kick velocity for the gas is «kick = (M*IMg)tc, 
or ~ 300 km s^^ for M, ~ Mg. There is thus certainly 
enough momentum input to make the gas highly turbu- 
lent. The cancellation of oppositely directed momentum re- 
sults in a reduction in this estimate. However, we have dis- 
cussed here only the momentum input; including supernova 
kinetic energy cou ld act to increase the above value (e.g., 
iDekel fc Silk. 19861 ). 

As we have shown in the paper, the effect of such feed- 
back is likely to result in some small angular momentum 
clouds or streams that may feed the SMBH. Our model 
therefore predicts a positive correlation of starburst activity 
with SMBH activity, with the latter being somewhat offset 
in time more often than not (it is also possible to activate 
an AGN earlier if there is some low angular momentum gas 
in the initial shell). 

Further detailed physics-based modeling of AGN feed- 
ing, feedback, and star formation, as well as comparisons to 
observations should shed light on the mode of quasar feed- 
ing, and also on the mode of bulge formation. 

8.4 Cosmological cold stream SMBH feeding? 

We have employed turbulent initial conditions as a simple 
and mathematically convenient way to introduce strong dis- 
ordered differential motions in the gas for controlled numeri- 
cal experiments in AGN feeding. We expect realistic galaxies 
to be always in the "turbulent" regime in the above sense be- 
cause the gas contracting onto a dark matter halo is unlikely 
to be of uniform density or possess a uniform velocity field. 
In fact, large-scale numerical simulations of galaxy forma- 
tion show cold gas streams that penetra te the shock-heate d 
media of massive dark matter haloes (jPekel et al.. 2009l l. 
This is physically similar to the ballistic motion of gas we ob- 



served in our simulations. While it is very unlikely that these 
very large-scale streams (hundreds of kpc) would be suffi- 
ciently well aimed to actually feed the SMBH, collisions of 
several of such streams in the centre of a galaxy would pump 
very strong differential motions and lead to some gas being 
set on small angular momentum orbits. We therefore spec- 
ulate that cold large-scale gas streams may promote AGN 
feeding by providing energy input to turbulence/differential 
motions on galactic scales. 



9 CONCLUSIONS 

The classical challenge to AGN feeding is the large an- 
gular momentum of gas expected to result in circu- 
larised discs of ~ kpc scales (e.g., Shlosman, Frank 
fc Begelman 1989). Transporting gas from these scales 
is extremely difficult due to long viscous times and 
the rapid consumption of gas in s tar formation (e.g., 
iGoodman. 200jNavakshin et al.. 20071 ). Here we pointed 
out that turbulence in the bulge, e.g., driven by supernova 
explosions, may be a way to overcome the AGN feeding de- 
bacle by setting some gas on ballistic orbits. 

To test these ideas, we have performed simulations of 
hydrodynamical gas flow with angular momentum, initially 
distributed in a spherical shell, in the static spherical poten- 
tial of a bulge with a central SMBH. The duration of our 
simulations is a few dynamical times of the shell, allowing us 
to study the formation of an accretion disc and the capture 
of gas within a small inner boundary region. For all values of 
the rotation velocity, the shell was not in equilibrium with 
the background potential. This was done deliberately, as the 
goal of this project was to study the formation of a disc from 
an infalling gas distribution from flrst principles. 

We found that angular momentum mixing in the runs 
without initial turbulence is very effictive, resulting in an ini- 
tial ring rather than an extended disc. Further evolution of 
the ring depends on whether it can transfer angular momen- 
tum quickly enough and avoid completely collapsing into 
stars. If star formation consumes the ring quickly then such 
rings are a dead end as far as SMBH feeding is concerned. 

Turbulent motions in the shell overcome the angular 
momentum mixing problem by creating high density regions 
that can travel nearly ballistically, retaining their initial an- 
gular momentum. Provided there is enough gas with small 
angular momentum the SMBH can be fed much more effi- 
ciently than in the "laminar" regime. The main conclusion 
here is that the star formation in the bulge may not simply 
deplete the gas (hence depriving the SMBH of its fuel) but 
may actually promote SMBH accretion by creating turbu- 
lent motions. 

We note that turbulence in this particular project is just 
a mathematically convenient formalism to introduce strong 
disordered differential motions in the gas. We expect that 
realistic galaxies are always in the turbulent regime because 
the state of the gas on large scales is unlikely to be that of 
uniform density and angular momentum. 

Further work with self-consistently driven turbulence 
from star formation and realistic initial (large-scale) bound- 
ary conditions is needed to establish the relevance of the 
"ballistic" accretion to feeding the SMBH, although our re- 
search suggests that it is a promising model. 
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APPENDIX A: AN ALTERNATIVE VELOCITY 
DISTRIBUTION FOR THE SHELL 

On the referee's request, we have re-simulated a fiducial pa- 
rameter set with an alternative velocity field. In this case we 
set I, the specific angular momentum, to be constant, rather 
than the rotation velocity. The magnitude of the velocity at 
a given radius is then 



where hot is a constant. This is a more conservative initial 
condition in the sense that this distribution would be ex- 
pected to form an infinitely thin ring, since the distribution 
is a delta function in angular momentum space. The broad- 
ening of this distribution due to turbulence would therefore 
be more difficult, and any enhancement to the accretion rate 
that this broadening provides is therefore a significant result. 
We note that the other "natural" choice of velocity field for 
a spherical shell, namely solid body rotation, is only reason- 
able for a system in hydrostatic equilibrium, which our shell 
is not. 

As can be seen from Figure lAll our main conclusion 
is unchanged for this alternative velocity distribution. The 
accretion rate onto the inner parsec is increased by ~ 2 or- 
ders of magnitude from zero/low turbulence to high turbu- 
lence. We note the the constant angular momentum within 
a spherical shell is, too, a strongly idealised setup, with sim- 
ilar behaviour at the polar regions to that of the constant 
Vrot condition. However, we are again attempting to model 
a simple case of gas mixing due to different angular momen- 
tum at similar radii, and unfortunately there is no "realistic" 
single-parameter velocity field that can describe a spherical 
distribution with these properties. 



APPENDIX B: CONVERGENCE TESTS ON 
NUMERICAL VISCOSITY AND RESOLUTION 

A great deal of the accretion that we find in our simu- 
lations is the result of turbulent motions creating dense, 
shocked gas. Due to the artificial viscosity present in the 
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Table 1. Initial conditions for each simulation. The ratio with respect to the gravitational potential energy gives an indication of how 
well supported the shell is against the external potential through virial motions: (-Eturb + ^thorm)/|E'grav | ~ 1/2 for a virialised shell. 
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*rcirc here refers to the initial value of circularisation radius for each of the particles, calculated in the region of the potential that they 
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Figure Al. Accretion rate versus time for a set of runs with 
a constant angular momentum velocity field, with Zrot = 0.02 
in code units. The other parameters for these runs were held at 
the fiducial values of S30, and color /linestyles is as standard for 
the other plots in the paper. The main result of the paper still 
holds: the accretion rate on the black hole strongly increases with 
increasing levels of turbulence when rotation is present (see il4.3|l . 

SPH code it is possible that some of this accretion is nu- 
merical rather than physical. We have therefore performed 
convergence tests both on the value of the artificial viscosity 
parameter, a (cf. Section [2]) and the numerical resolution of 
our fiducial simulations S30 (no turbulence, Viot = 0.3) and 
S35 (i)turb = 1, "Viot — 0.3), comparing the accretion rate 
between the two. The results are shown in Figures IBll and 
IB2l for the artificial viscosity and particle resolution respec- 
tively. 

For Figure IBll there is clear convergence both in the 
runs with no turbulence and those with «turb = 1. This 
suggests that the role of the artificial viscosity is negligible in 
determining the accretion rate through the SMBH accretion 
radius, at least until the main ballistic mode of feeding ends 
(it is likely that once the accretion becomes dominated by 
an accretion disc mode, i.e., to late times, the numerical 
viscosity will play a greater role). 

Figure [B2l shows the convergence results on the number 
of SPH particles. In the runs with no turbulence it is not 
clear whether the accretion rate has converged with increas- 
ing resolution, as we were unable to run simulations with a 
yet higher number of particles due to computational limita- 
tions. However, the plot would seem to suggest that by the 
highest resolution there is some convergence to late times. 
It may of course simply be the case that for this particular 
value of Wrot all of the accretion is numerical for the runs 
with zero turbulence, and the accretion rate will therefore 
vanish in the limit of increasing resolution. We note that 
this leaves our conclusions unaffected, and indeed strength- 
ens the case for turbulence enhancing the accretion rate onto 
the SMBH. In the runs with uturb = 1 there is some reason- 
able convergence, albeit about a relatively large spread with 
no clear trend. A possible reason for this is the fact that the 
SPH method tends to suffer from an inability to converge 
with increasing resolution in regions where particles are dis- 



Figure Bl. Accretion rate versus time for convergence tests on 
the value of the artificial viscosity parameter, a. Solid linestyles 
refer to simulations without turbulence while dashed linestyles 
refer to simulations with t>turb = 1- The rotation velocity for all 
runs here was iirot = 0.3. Colors are: a = 0.1 (green), a = 0.3 
(blue), a = 0.5 (red), a. = 0.7 (magenta) and a = 1.0 (black; our 
fiducial value for the simulations in the paper). The simulations 
with no seeded turbulence show good convergence with a trend 
towards an increased accretion rate with an increasing a, while 
the simulations with I'turb = 1 show good convergence but with 
no clear trend. 
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Figure B2. Accretion rate versus time for convergence tests on 
the numerical resolution. Once again solid linestyles are for runs 
with i)turb = while dashed refer to runs with Dturb = Ij a-nd 
all simulations shown used Vrot = 0.3. Colors refer to A^sph = 
5.4 X 10^ (green), AfspH = 8.0 X 10^ (blue), AfspH = 1.2 X 10'^ 
(magenta), A^sfh = 1-8 x 10"^ (red), A^SFH = 2.7 X 10^ (brown) 
and Afgpjj = 4 X 10® (black; our fiducial resolution). In the runs 
with no seeded turbulence there is a strong trend of decreasing 
accretion with increasing resolution, which shows signs of starting 
to converge by the time the highest resolution is reached. In the 
runs with "Uturb = 1 there is convergence to within a factor of a 
few but without a clear trend in the accretion rate. 
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tributed irregularly on t he kernel scale (su ch as in strong 
shocks), as discussed by iRead et al. However, this 

spread is not large enough to cast any doubt on our main 
conclusion, and furthermore is only ^ a factor of a few for 
the 3 higher resolution runs (red, brown, black lines). 

We can therefore be confident that the effect seen in 
our simulations of an enhanced accretion rate with increas- 
ing turbulence (up to saturation) through a ballistic feeding 
mode is physical rather than numerical. 



